A máquina de escrever do meu avô
…da maldição do erro operacional
Fonte: (QLD 2016)
PS.: ele não é Nobel…
…e ainda…
Fonte: (Boddy 2016)
Meta final: que vários dos processos de análise de dados da empresa passem a ser feitos dentro do fluxo de trabalho do R.
Fonte:(Frigaard 2017)
Ao fim do curso o objetivo é que todos os fios da meada sejam puxados para que o aluno consiga continuar por si só usando a vasta documentação disponível.
Fonte:(Wickham and Grolemund 2016)
Fonte: (Frigaard 2017)
(Mello 2019)
R é uma linguagem que é interpretada por um engine gratuito.
RStudio é o melhor ambiente de programação da linguagem R. A versão mais simples, que é totalmente funcional, é gratuita.
Na visualização padrão, ele oferece um console para execução de comandos e uma janela com a visualização dos environments, ou seja, das variáveis que ele guarda na sessão atual.
Existe orientação a objetos no R, mas não está no escopo deste curso
Note que não há variáveis que armazenam dado escalar, como já vimos.
Dentre os vetores há:
vetores atômicos (seus elementos são do mesmo tipo primário)
listas (seus elementos, que são vetores atômicos, são de tipos primários diferentes)
Tipos de vetores
Fonte: (Wickham 2019)
Os vetores atômicos podem ser dos seguintes tipos:
Tipos primários
Fonte: (Wickham 2019)
A biblioteca bench oferece funções ótimas para avaliar a performance de pedaços pequenos de código.
resultados_perf <- mark(
sem_loop(1e4),
com_loop(1e4),
(1:1e4)^2
)
#aprenderemos o que é %>% e select() posteriormente
resultados_perf %>%
select(expression, min, median, `itr/sec` )## # A tibble: 3 x 4
## expression min median `itr/sec`
## <bch:expr> <bch:tm> <bch:tm> <dbl>
## 1 sem_loop(10000) 6.91ms 7.41ms 127.
## 2 com_loop(10000) 129.33ms 140.05ms 7.17
## 3 (1:10000)^2 17.43us 19.53us 18415.
Monty Hall era uma espécie de Sílvio Santos juvenil (sub 80) americano.
Um dos seus jogos consistia em mostrar três portas ao otár… (ops) convidado. Em uma delas tem um carro.
Antes do resultado, o apresentador revela uma das portas e pergunta se o convidado que trocar a escolha.
O que vocês acham? Melhor trocar, manter a escolha original ou tanto faz?
Vamos ver… mas dá pra simular sem saber quase nada.
Vamos usar uma das funções da família r<familia de distribuição de prob>(). Neste caso, a rbinom, que simula a distribução binomial (aquela que equivale ao evento de jogar n moedas (ou alguma coisa com dois lados) para cima e ver quantas deram cara).
n_simul <- 10000
n_questoes <- 240
min_aprovacao <- 0.6
n_aprovado <- 240 * min_aprovacao
prob_questao <- 0.2
acertos <- rbinom(n = n_simul, size = n_questoes, prob = prob_questao )
sum(acertos >= n_aprovado)/n_simul ## [1] 0A chance é praticamente nula.
Na verdade, a grande massa da distribuição fica muito distante.
dado <- enframe(acertos/n_questoes)
mostra_chances <- function(acertos, n_questoes){
ggplot(enframe(acertos/n_questoes)) +
geom_density( aes(x = value)) +
scale_x_continuous(
labels = percent_format(accuracy = 1),
limits = c(0,1),
breaks = seq(0, 1, 0.1)
) +
labs(x ="% Acertos") +
geom_vline(xintercept = min_aprovacao, color = "red") +
theme_light()
}
mostra_chances(acertos, n_questoes)
O exemplo anterior era muito simplista: ninguém chuta tudo.
Imagine que sabemos qual a chance de aparecer uma pergunta onde podemos descartar 0 alternativas, a chance de uma onde descartamos 1 e assim por diante.
#definindo a chance podermos eliminar 0, 1, 2, ... 4 alternativas
fracao_eliminar_questoes <- c( 0.1, 0.1, 0.2, 0.25 , 0.35 )
#definindo o número de questões
n_questoes_cada_elimina <- t(rmultinom(n_simul, size = n_questoes, fracao_eliminar_questoes))
probs_quando_elimina <- 1/(5:1)
acertos_concatenados <-
rbinom(
n = n_simul * 5 ,
size = as.vector(t(n_questoes_cada_elimina)),
prob = probs_quando_elimina
)## [,1] [,2] [,3] [,4] [,5]
## [1,] 23 24 50 47 96
## [2,] 24 17 48 64 87
## [3,] 23 30 43 56 88
## [4,] 20 25 51 52 92## [1] 4 6 17 25 96 4 3 13 34 87 6 5 13 36 88 4 9 17 25 92## [,1] [,2] [,3] [,4] [,5]
## [1,] 4 6 17 25 96
## [2,] 4 3 13 34 87
## [3,] 6 5 13 36 88
## [4,] 4 9 17 25 92
## [5,] 6 5 11 30 92## [1] 0.3131
Arrumar os dados de forma que as linhas sejam eventos e as colunas sejam atributos do evento ajuda muito a rodar modelos e construir visualizações eficientemente.
O que é o evento e o que é o atributo pode variar até para diferentes usos do mesmo dado. Mas a prática ajuda a determinar isso.
dplyr
Fonte: (Courtiol 2019)
rename()rename() é usada para modificar os nomes das colunas. Ela renomeia as colunas indicadas e mantném as outras.
A função mutate é usada para criar novas colunas no tibble, mantendo as outras. 
mutate() (cont.):Outras funções úteis são as que fazem operações acumuladas e as operações de lag() e lead()
## Registered S3 method overwritten by 'xts':
## method from
## as.zoo.xts zoo## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoo## Registered S3 methods overwritten by 'forecast':
## method from
## fitted.fracdiff fracdiff
## residuals.fracdiff fracdiff##
## BETS-package: Found 33 out of 18706 time series.## # A tibble: 33 x 7
## code description unit periodicity start last_value source
## <chr> <chr> <chr> <chr> <chr> <chr> <chr>
## 1 1 Exchange rate - Free~ c.m.u~ D 28/11~ 26/03/2018 Sisbac~
## 2 10813 Exchange rate - Free~ c.m.u~ D 28/11~ 03/05/2018 Sisbac~
## 3 11753 Real effective excha~ Index M 31/01~ mar/2018 BCB-DS~
## 4 11758 Real effective excha~ Index M 31/01~ mar/2018 BCB-DS~
## 5 11763 Real effective excha~ Index M 31/01~ mar/2018 BCB-DS~
## 6 11768 Real effective excha~ Index M 31/01~ mar/2018 BCB-DS~
## 7 17887 Exchange rate (perio~ US$/c~ M 31/01~ nov/2016 IMF-IFS
## 8 17888 Exchange rate (perio~ US$/c~ M 31/01~ nov/2016 IMF-IFS
## 9 18441 Exchange rate (end o~ US$/c~ M 31/01~ nov/2016 IMF-IFS
## 10 18442 Exchange rate (end o~ US$/c~ M 31/01~ nov/2016 IMF-IFS
## # ... with 23 more rowsNo código abaixo, calculamos o retorno da série, a volatilidade histórica e a volatilidade EWMA
## Registered S3 method overwritten by 'xts':
## method from
## as.zoo.xts zoo## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoo## Registered S3 methods overwritten by 'forecast':
## method from
## fitted.fracdiff fracdiff
## residuals.fracdiff fracdiffdolar_com_vol <- dolar %>%
filter(date > ymd("1994-07-01")) %>%
arrange(date) %>%
mutate(
retorno = (value - lag(value))/value,
retorno_quad = retorno^2,
dia = row_number(),
fator_ewma = (1/0.94)^dia*1e-20
) %>%
filter(!is.na(retorno)) %>%
mutate(vol = sqrt(cumvar(retorno)) * sqrt(252) ) %>%
mutate(vol_ewma = sqrt(cumsum(retorno_quad * fator_ewma)/cumsum(fator_ewma)) * sqrt(252) ) %>%
rename(dolar = value) %>%
select(
date,
dolar,
retorno,
vol,
vol_ewma
)
datatable(dolar_com_vol) %>%
formatPercentage(c("retorno", "vol", "vol_ewma"), 2)dolar_ajeitado <- dolar_com_vol %>%
gather(variavel, valor, - date)
dolar_ajeitado %>%
ggplot() +
geom_line(aes(x = date, y = valor)) +
facet_grid( variavel ~ . , scales = "free") +
theme_light() 
transmute():transmute() cria colunas e mantém apenas as colunas citadas
## # A tibble: 6 x 3
## ano pais pib
## <int> <fct> <dbl>
## 1 1952 Afghanistan 6567086330.
## 2 1957 Afghanistan 7585448670.
## 3 1962 Afghanistan 8758855797.
## 4 1967 Afghanistan 9648014150.
## 5 1972 Afghanistan 9678553274.
## 6 1977 Afghanistan 11697659231.A função arrange serve para ordenar as linhas do tibble.
filter():filter() seleciona colunas de acordo com os seus valores
filter() (cont.):Filtrando países (note o operador %in%)
gapminder %>%
filter(country %in% c("Brazil", "Argentina", "Chile")) %>%
ggplot() +
geom_line(aes(x = year, y = gdpPercap, color = country )) +
geom_point(aes(x = year, y = gdpPercap, color = country )) +
theme_light()
top_n():top_n seleciona as n linhas maiores de acordo com uma das colunas.
top_n() (cont.):Selecionando os países mais ricos em 2007.
Depois aprenderemos como ordenar essas barras
gapminder %>%
filter(year == 2007) %>%
top_n(5, gdpPercap) %>%
ggplot() +
geom_col(aes(x = country, y = gdpPercap)) +
theme_light() 
top_frac():
top_frac():ricos <- gapminder %>%
filter(year == 2007) %>%
top_frac(.2, gdpPercap ) %>%
mutate(categoria = "Rico")
pobres <- gapminder %>%
filter(year == 2007) %>%
top_frac(.2, desc(gdpPercap) ) %>%
mutate(categoria = "Pobre")
bind_rows(ricos, pobres) %>%
ggplot(aes(x = gdpPercap, y = lifeExp)) +
geom_point( aes(color = continent)) +
facet_grid(. ~ categoria, scales = "free_x") +
geom_smooth(method = "lm", se = FALSE) +
theme_light()
slice():
group_by():
summarise() (cont.):
group_by() (cont.):A função group_by só leva a uma sumarização, ou seja, só transforma o tibble em um tibble com o número de linhas igual ao número de grupos, quando executamos a função summarise()
Se executarmos summarise sem particionar o tibble, a operação resulta em uma linha.
maiores_temp_dia <- dados_heathrow %>%
group_by(date(date)) %>%
summarise(
maxima = max(air_temp),
minima = min(air_temp),
media = mean(air_temp)
)
datatable(maiores_temp_dia) %>%
formatRound(c("maxima", "minima", "media"), 1)A função top_n retorna os n maiores valores. Se o tibble estiver agrupado, pra cada grupo.
maiores_temp_dia <- dados_heathrow %>%
group_by(date(date)) %>%
top_n(1, air_temp) %>%
ungroup() %>%
mutate(
hora = hour(date),
estacao =
case_when(
month(date) %in% 1:3 ~ "Inverno",
month(date) %in% 7:9 ~ "Verão",
TRUE ~ "Outono/Primavera"
)
) %>%
select(hora, estacao, air_temp)
ggplot(maiores_temp_dia) +
geom_density( aes(x = hora, color = estacao )) +
theme_light()
O caso mais comum é ler dados em formato de tabela para um tibble
Fonte: (RStudio 2019a)
Uma página WEB pode ser representada por uma árvore de objetos, também chamada de DOM (Document Object Model).
Esta árvore de objetos é definida pelo conteúdo de linguagem html que existe na página (e pode ser modificado por scripts em javascript e definições de estilo do CSS).
Podem existir objetos de vários tipos em uma página: links, inputs de dados, tabelas, células de tabelas, parágrafos, cabeçalhos etc.
Para retirar da página web o conteúdo de que precisamos, temos que analisar como é esta árvore de objetos e que nós desta árvore nos interessam.
Imagine que queremos buscar dados na página de histórico de preços e taxas dos títulos brasileiros
A tecla F12 do Chrome nos permite ver a árvore DOM da página em que estamos navegando.
É possível clicar com o botão direito e inspecionar um elementos específico de forma a saber onde ele está na árvore e que tipo de elemento ele é (mesmo que você saiba pouco de html).
O mais importante é saber que uma tag html que define um elemento tem a sintaxe:
<tipo_elemento nome_atributo=valor_attributo>texto do elemento</tipo_elemento>
Mesmo sem saber hmtl, fica claro que queremos esse tal de atributo href dos tais elementos a seja lá o que diabos isso seja (a é um link e href é o destino do link).
Aí fica fácil fazer a análise que desejarmos
ntnb_2045 <- taxas_titulos %>%
filter(
titulo == "NTN-B Princ",
vencimento == ymd("2045-05-15")
) %>%
mutate(taxa = (taxa_compra_manha +taxa_venda_manha)/2 )
ggplot(ntnb_2045) +
geom_line(aes(x = data, y = taxa) ) +
scale_x_date(
date_breaks = "3 months",
limits = c(ymd("2016-12-01"),NA),
labels = date_format("%m/%y")
) +
scale_y_continuous(
labels = percent_format(accuracy = 0.1)
) +
labs(y = "NTN-B Principal 2045", x = "Data") +
theme_light()
Nem todas as páginas são fáceis de ler.
A página que mostra os DIs na BMF, por exemplo é esquisita:
O nosso objeto de interesse aparece na ferramenta chamada por F12, mas não é encontrada pela rvest ao ler o HTML
Existem duas técnicas para o Web Scraping:
Simular a navegação num browser emulado (fora do escopo deste curso)
Interpretar páginas com base no conteúdo recebido
Na página que visualizamos da BMF, o conteúdo recebido inclui um script javascript que popula a tabela, por isso não conseguimos reconhecer o conteúdo da tabela de pronto, pois ela só é carregada pela execução do script.
Conseguimos, no entanto, extrair do script as informações de que precisamos
Reparamos que o data frame no slide anterior não está “Tidy”, ou seja não está de forma que cada linha represente um evento e cada coluna represente um atributo do evento.
Para nós, neste caso, um evento é formado por todas as informações de um ativo em um dia e não uma só das informações de um ativo em um dia.
Isso porque é extremamente comum fazermos contas com mais de uma informação do ativo em um dia (máxima - mínima, por exemplo)
O pacote Tidyr ajuda a arrumar os data frames dessas formas. O hex sticker dele é bem explicativo.
pivot_longer() e pivot_wider()Os nomes das principais funções mudaram em setembro de 2019 (quando saiu a versão 1.0.0). Antes se chamavam gather() e spread() e agora se chamam pivot_wider() e pivot_longer(), o que é mais intuitivo.
Para juntar as informações de dois tibbles em um só, podemos fazer isso de três formas
Fonte: (RStudio 2019b)
Fonte: (Exame 2018)
Às vezes mesmo sem querer (será?)
PS. gráfico de 2014
Fonte: (Alves 2014)
Gráfico com sombra, 3D, estilos que dificultam a interpretação
Fonte: (Healy 2018)
O pessoal que usa Excel muitas vezes ama 3D, mas…
(Healy 2018)
Fonte: (White 2015)
Assim como os restaurantes devemos servir pizzas de dois sabores no máximo
Tentem descobrir qual a menor e a menor categoria em cada caso
Fonte: (Viz 2018)
Fonte: (Viz 2018)
Gráficos de dois eixos podem mostrar relações espúrias muito facilmente. E elas dependem da escolha da escala e dos limites dos eixos.
Fonte (Rost 2018)
O mesmo dado pode levar aos seguintes gráficos (e suas soluções)
Enfatize o dado relevante
Integre texto e gráfico
Use fontes e cores diferentes do padrão
Principalmente em colunas e barras, faça o eixo começar de ZERO
loadfonts(device = "win")
by_country <- organdata %>% group_by(consent_law, country) %>%
summarize_if(is.numeric, list(mean = mean, sd = sd), na.rm = TRUE) %>%
ungroup()
p <- ggplot(data = by_country,
mapping = aes(x = gdp_mean, y = health_mean))
p + geom_point( color = "coral4") +
geom_text_repel(data = subset(by_country, gdp_mean > 25000),
mapping = aes(label = country),
size = 3,
color = "coral4",
family="Bookman Old Style"
) +
labs(y = "Gastos com saúde per cap.", x = "PIB per cap." ) +
scale_x_continuous(
labels = number_format(decimal.mark = ",", big.mark = "."),
limits = c(0,NA)
) +
scale_y_continuous(
labels = number_format(decimal.mark = ",", big.mark = "."),
limits = c(0,NA)
) +
theme_minimal(
) +
ggtitle(
label = "Gastos em saúde x PIB"
) +
theme(
text=element_text(family="Bookman Old Style", color = "coral4"),
axis.text = element_text(colour = "coral4"),
panel.background = element_rect(fill = "beige"),
panel.grid.minor = element_line(colour = "bisque1"),
panel.grid.major = element_line(colour = "bisque3")
)
Use small multiples: divida o gráfico em gráficos iguais menores cada um com parte dos dados, seguindo uma regra categórica
#("azure2", "chocolate", "steelblue4", "darkslategray", "slategray3", "slategray4")
p <- ggplot(data = gapminder, mapping = aes(x = year, y = gdpPercap))
p + geom_line(color="chocolate", aes(group = country)) +
geom_smooth(size = 1.1, method = "loess", se = FALSE, color = "tan4") +
scale_y_log10(labels=scales::dollar) +
facet_wrap(~ continent, ncol = 5) +
labs(x = "Year",
y = "GDP per capita",
title = "GDP per capita on Five Continents")+
theme_minimal() +
theme(
text=element_text(family="CMU Serif", color = "darkslategray", size = 8),
axis.text = element_text(colour = "darkslategray"),
panel.background = element_rect(fill = "azure2"),
panel.grid.minor = element_line(colour = "slategray2"),
panel.grid.major = element_line(colour = "slategray2"),
strip.text = element_text(family="CMU Serif", colour = "darkslategray"),
aspect.ratio = 1
)
A GGPLOT2 parece mais complicada de usar do que funções como a plot() da biblioteca base do R.
Talvez porque as pessoas não são sempre apresentadas às camadas da gramática “Grammar of Graphics” que fundamenta a biblioteca.
Fonte: (Scavetta 2018)
Fonte: (Scavetta 2018)
Veja que, pela variável ser representada de forma discreta (uma casa decimal), há sobreposição de pontos. para dar a real impressão de quantos pontos existem, o ideal é inserir um ruído.
ggplot(iris, aes(x = Sepal.Length, y = Sepal.Width)) +
geom_jitter(alpha = 0.6) +
facet_grid(. ~ Species)
ggplot(iris, aes(x = Sepal.Length, y = Sepal.Width)) +
geom_jitter(alpha = 0.6) +
facet_grid(. ~ Species) +
stat_smooth(method = "lm", se = F, col = "red")
ggplot(iris, aes(x = Sepal.Length, y = Sepal.Width)) +
geom_jitter(alpha = 0.6) +
facet_grid(. ~ Species) +
stat_smooth(method = "lm", se = F, col = "red") +
coord_flip()
wes <- wes_palette("Royal2", 5, "discrete")
ggplot(iris, aes(x = Sepal.Length, y = Sepal.Width)) +
geom_jitter(alpha = 0.6, color = wes[1]) +
facet_grid(. ~ Species) +
stat_smooth(method = "lm", se = F, col = wes[1]) +
coord_flip() +
theme_minimal() +
theme(
text=element_text(family="Bradley Hand ITC", color = wes[1], size = 16),
axis.text = element_text(colour = wes[1]),
panel.background = element_rect(fill = "white"),
panel.grid.minor = element_line(colour = wes[2]),
panel.grid.major = element_line(colour = wes[3]),
strip.text = element_text(family="Bradley Hand ITC", color = wes[1]),
panel.border = element_rect(colour = wes[4], fill = NA)
)
gols <- read_csv("dados/football_events/ginf.csv") %>%
select(fthg, ftag) %>%
pivot_longer(cols = everything(), names_to = "casa_fora", values_to = "gols")
ggplot(gols) +
geom_histogram(
aes(x = gols),
binwidth = 1
) +
scale_x_continuous(breaks = 0:10, minor_breaks = c()) +
theme_minimal()
result <- fitdistr( gols$gols , densfun = "Poisson" )
pois <- rpois(length(gols$gols ), lambda = result$estimate ) %>%
enframe(value = "gols") %>%
mutate(tipo = "poisson")
real_simulado <- gols %>%
mutate(tipo = "real" ) %>%
bind_rows(pois)
ggplot() +
geom_histogram(
data = real_simulado,
aes(x = gols, group = tipo, fill = tipo),
binwidth = 1,
show.legend = TRUE,
position = "identity",
alpha = 0.5
) +
scale_x_continuous(breaks = 0:10, minor_breaks = c()) +
geom_vline(aes(xintercept = result$estimate)) +
geom_text(aes(x = result$estimate, y = 7000 ),nudge_x = 0.3, label = number(result$estimate, accuracy = 0.01 )) +
theme_minimal()
O gráfico mostrando a função de probablidade acumulada empírica mostra propriedades que o histograma não mostra
ggplot(gols, aes(x = gols)) +
stat_ecdf(geom = "step") +
scale_x_continuous(breaks = 0:10, minor_breaks = c()) +
theme_minimal()
Para variáveis contínuas, faz sentido usar o gráfico de densidade
ufc_pesos <- read_csv("dados/ufc/data.csv") %>%
filter(weight_class == "Heavyweight") %>%
select(
R_Height_cms,
B_Height_cms,
Winner
) %>%
transmute(
altura_vencedor = if_else(Winner == "Red", R_Height_cms, B_Height_cms ),
altura_perdedor = if_else(Winner == "Red", B_Height_cms, R_Height_cms )
) %>%
pivot_longer(cols = everything(), names_to = "resultado", values_to = "altura")
ggplot(ufc_pesos) +
geom_density(aes(x = altura)) +
theme_minimal()
Mais uma vez a distribuição empírica nos deixa ver mais coisa: os quantis
## Warning: Removed 1 rows containing non-finite values (stat_ecdf).
## Warning: Removed 1 rows containing non-finite values (stat_density).
ggplot(ufc_pesos) +
stat_ecdf(aes(x = altura, color = resultado), geom = "density" ) +
theme_minimal()## Warning: Removed 1 rows containing non-finite values (stat_ecdf).
ufc_pesos <- read_csv("dados/ufc/data.csv") %>%
select(
weight_class,
R_Height_cms,
B_Height_cms,
Winner
) %>%
transmute(
weight_class,
altura_vencedor = if_else(Winner == "Red", R_Height_cms, B_Height_cms ),
altura_perdedor = if_else(Winner == "Red", B_Height_cms, R_Height_cms )
) %>%
pivot_longer(cols = -weight_class, names_to = "resultado", values_to = "altura")
ggplot(ufc_pesos) +
geom_density_ridges(aes(x = altura, y = weight_class ), fill = "lightblue") +
theme_ridges() +
scale_x_continuous(breaks = seq(150,by = 5, to = 220))
Legal a ideia, mas não queremos comparar mulheres também
A biblioteca stringrfacilita muio lidar com strings. Por exemplo detectar um pedaço de string em outra
ufc_pesos <- read_csv("dados/ufc/data.csv") %>%
filter(!str_detect(weight_class, "Women")) %>%
select(
weight_class,
R_Height_cms,
B_Height_cms,
Winner
) %>%
transmute(
weight_class,
altura_vencedor = if_else(Winner == "Red", R_Height_cms, B_Height_cms ),
altura_perdedor = if_else(Winner == "Red", B_Height_cms, R_Height_cms )
) %>%
pivot_longer(cols = -weight_class, names_to = "resultado", values_to = "altura")
ggplot(ufc_pesos) +
geom_density_ridges(aes(x = altura, y = weight_class ), fill = "lightblue") +
theme_ridges() +
scale_x_continuous(breaks = seq(150,by = 5, to = 220))
Melhorou, mas ainda á estranho. Melhor seria se as classes de peso estivessem em ordem
No R, as variáveis categóricas são chamadas de fatores. Cada valor possível de um fator é representadas internamente com um identificador numérico, e não com a própria string.
Cada valor possível da variável categórica também é chamada de level
Podemos manipular essas representações facilmente com a forcats.
No caso, queremos ordenar nosso fator com a ordem as classes de peso no UFC
A função fct_relevel possibilita a ordenação manual de uma variável factor
ufc_pesos_ordem <- ufc_pesos %>%
mutate(
weight_class = fct_relevel(weight_class,
c(
"Flyweight",
"Bantamweight",
"Featherweight",
"Lightweight",
"Welterweight",
"Middleweight",
"Light Heavyweight",
"Heavyweight",
"Catch Weight",
"Open Weight"
)
)
)
ggplot(ufc_pesos_ordem) +
geom_density_ridges(aes(x = altura, y = weight_class ), fill = "lightblue") +
theme_ridges() +
scale_x_continuous(breaks = seq(150,by = 5, to = 220))
Os gráficos do tipo scatter plot
Vamos brincar um pouco com os dados do Banco Mundial
taxa_homicidios <- wb("VC.IHR.PSRC.P5", country = "all", mrv = 10, POSIXct = TRUE ) %>%
select(iso3c, value, date) %>%
group_by(iso3c) %>%
top_n(1, date) %>%
rename(homicidios = value)
gini <- wb("SI.POV.GINI", country = "all", mrv = 10, POSIXct = TRUE ) %>%
select(iso3c, value, date) %>%
group_by(iso3c) %>%
top_n(1, date) %>%
rename(desigualdade = value)
pib_per_capita <- wb("NY.GDP.PCAP.PP.KD", country = "all", mrv = 10, POSIXct = TRUE ) %>%
select(iso3c, value, date) %>%
group_by(iso3c) %>%
top_n(1, date) %>%
rename(riqueza = value)
pib_gini_homi <- taxa_homicidios %>%
inner_join(gini, by = c("iso3c")) %>%
inner_join(pib_per_capita, by = c("iso3c")) %>%
left_join(codelist, by = c("iso3c") ) %>%
select(riqueza, desigualdade, homicidios, continent, region) %>%
pivot_longer(cols = c("riqueza", "desigualdade"), names_to = "tipo_indice", values_to = "indice" ) %>%
mutate(tipo_indice = str_to_title(tipo_indice))ggplot(pib_gini_homi, aes(y = homicidios, x = indice)) +
geom_point(alpha = 0.4) +
facet_wrap(~tipo_indice, scales = "free", ncol = 1) +
geom_smooth() +
theme_minimal() +
labs(x = "", y = "Homicídios por 100 mil") +
scale_x_continuous(labels = number_format(decimal.mark = ",", big.mark = ".")) 
É possível visualizar os eixos em log na base 10
ggplot(pib_gini_homi, aes(y = homicidios, x = indice)) +
geom_point(alpha = 0.4) +
facet_wrap(~tipo_indice, scales = "free", ncol = 1) +
geom_smooth() +
theme_minimal() +
labs(x = "", y = "Homicídios por 100 mil") +
scale_x_continuous(labels = number_format(decimal.mark = ",", big.mark = ".")) +
scale_y_log10() 
É possível também usar as cores para olhar uma terceira característica
ggplot(pib_gini_homi, aes(y = homicidios, x = indice)) +
geom_point(alpha = 0.4, aes(color = continent )) +
facet_wrap(~tipo_indice, scales = "free", ncol = 1) +
geom_smooth() +
theme_minimal() +
labs(x = "", y = "Homicídios por 100 mil") +
scale_x_continuous(labels = number_format(decimal.mark = ",", big.mark = ".")) +
scale_y_log10() +
theme(legend.position = "top")
A linha de regressão também pode ser estimada para cada grupo. Neste caso vamos usar o método de regressão linear, pois o LOESS original nao funciona nada bem com muito poucos dados (nem a linear, mas…)
ggplot(pib_gini_homi, aes(y = homicidios, x = indice, color = continent)) +
geom_point(alpha = 0.4 ) +
facet_wrap(~tipo_indice, scales = "free", ncol = 1) +
geom_smooth(se = FALSE, method = "lm") +
theme_minimal() +
labs(x = "", y = "Homicídios por 100 mil") +
scale_x_continuous(labels = number_format(decimal.mark = ",", big.mark = ".")) +
scale_y_log10() 
analise_desigualdade <- pib_gini_homi %>%
filter(tipo_indice == "Desigualdade")
ggplot(analise_desigualdade, aes(y = homicidios, x = indice, color = continent)) +
geom_point(alpha = 0.4) +
facet_wrap(~tipo_indice, ncol = 1) +
geom_smooth(se = FALSE, method = "lm") +
theme_minimal() +
labs(x = "Índice de Gini", y = "Homicídios por 100 mil") +
scale_x_continuous(labels = number_format(decimal.mark = ",", big.mark = ".")) +
scale_y_log10() +
facet_wrap(~continent) +
theme(legend.position = "top")
ggplot(analise_desigualdade, aes(y = homicidios, x = indice, color = continent)) +
geom_text(aes(label = iso3c), size = 2) +
facet_wrap(~tipo_indice, ncol = 1) +
geom_smooth(se = FALSE, method = "lm") +
theme_minimal() +
labs(x = "Índice de Gini", y = "Homicídios por 100 mil") +
scale_x_continuous(labels = number_format(decimal.mark = ",", big.mark = ".")) +
scale_y_log10() +
facet_wrap(~continent) +
theme(legend.position = "top")
Podemos avaliar três variáveis contínuas usando um gráfico de bolhas
Primeiro vamos colocar em wide as duas variáveis que tínhamos deixado long.
pib_gini_homi_wide <- pib_gini_homi %>%
pivot_wider(names_from = tipo_indice, values_from = indice )
head(pib_gini_homi_wide)## # A tibble: 6 x 6
## # Groups: iso3c [6]
## iso3c homicidios continent region Riqueza Desigualdade
## <chr> <dbl> <chr> <chr> <dbl> <dbl>
## 1 ALB 2.3 Europe Southern Europe 12306. 29
## 2 DZA 1.4 Africa Northern Africa 13886. 27.6
## 3 AGO 4.8 Africa Middle Africa 5725. 42.7
## 4 ARG 5.1 Americas South America 18282. 41.2
## 5 ARM 2.4 Asia Western Asia 9178. 33.6
## 6 AUS 0.8 Oceania Australia and New Zealand 45439. 35.8A escala de cor usada abaixo, Viridis, oferece a mesma sensação para olorido e preto e branco e é feita para ajudar daltônicos.
O log foi usado aqui para evitar que os valores extremos dificultem a visualização da escala de cores.
ggplot(pib_gini_homi_wide) +
geom_point(aes(color = homicidios, x = Desigualdade , y = Riqueza )) +
scale_color_viridis_c(trans = "log10", direction = -1, option = "inferno") +
theme_minimal()
hdi <- read_csv("dados/hdi/atlas.csv") %>%
rename(municipio = `município` ) %>%
select(
ano,
rdpc,
gini,
municipio,
uf
) %>%
filter(ano == 2010)
violencia <- extract_text("dados/atlasviolencia/8099-tabelamunicipiostodossite.pdf", encoding = "UTF-8") %>%
str_split(pattern = fixed("\n")) %>%
unlist() %>%
enframe( value = "linha") %>%
mutate(
UF = str_extract(linha, "[A-Z]{2}") ,
municipio = str_extract(linha, "(?<=[A-Z]{2} ).+?(?=[0-9])"),
numeros = str_extract_all(linha, "(?<= )[0-9 ,.]*")
) %>%
unnest_legacy() %>%
filter(str_length(numeros)>1) %>%
mutate(
numeros = str_trim(numeros),
municipio = str_trim(municipio) %>% str_to_upper()
) %>%
separate(
col = numeros,
into = c("populacao", "homicidios", "ocultos", "taxa_homicidios"),
sep = " "
) %>%
mutate(taxa_homicidios = parse_number(taxa_homicidios, locale = locale(decimal_mark = ",")) )
ufs <- municipios_ajeitado %>%
select(
UF,
cod_municipio
) %>%
mutate(
cod_uf = cod_municipio %/% 100000
) %>%
select(-cod_municipio) %>%
distinct() %>%
filter(!is.na(cod_uf))
hdi %<>% inner_join(ufs, by = c("uf" = "cod_uf"), suffix = c("_old", ""))
hdi_violencia <- hdi %>%
inner_join(violencia, by = c("municipio", "UF") ) %>%
mutate(
regiao = case_when(
uf %/% 10 == 1 ~ "Norte",
uf %/% 10 == 2 ~ "Nordeste",
uf %/% 10 == 3 ~ "Sudeste",
uf %/% 10 == 4 ~ "Sul",
uf %/% 10 == 5 ~ "Centro-Oeste",
TRUE ~ NA_character_
)
)ggplot(hdi_violencia) +
geom_jitter(aes(color = taxa_homicidios, x = gini , y = rdpc ), alpha = 0.2 ) +
scale_color_gradient(low = "blue", high = "red", trans = "log10" ) +
facet_wrap(~regiao) +
theme_minimal()## Warning: Transformation introduced infinite values in discrete y-axis
ufc_data <- read_csv("dados/ufc/data.csv") %>%
select(weight_class, no_of_rounds )
ufc_raw_data <- read_csv2("dados/ufc/raw_total_fight_data.csv") %>%
select(last_round)
ufc_tudo <- bind_cols(ufc_data, ufc_raw_data) %>%
filter(no_of_rounds == 3) %>%
group_by(weight_class) %>%
mutate(n_classe = n()) %>%
group_by(weight_class, last_round) %>%
summarise(n_round = n(), n_classe = mean(n_classe)) %>%
mutate(frac_lutas = n_round/n_classe ) %>%
ungroup() %>%
filter( weight_class %in%
c(
"Flyweight",
"Bantamweight",
"Featherweight",
"Lightweight",
"Welterweight",
"Middleweight",
"Light Heavyweight",
"Heavyweight"
)) %>%
mutate(
weight_class = fct_relevel(weight_class,
c(
"Flyweight",
"Bantamweight",
"Featherweight",
"Lightweight",
"Welterweight",
"Middleweight",
"Light Heavyweight",
"Heavyweight"
)
)
)ggplot(ufc_tudo, aes(x = last_round, y = weight_class)) +
geom_tile(aes(fill = frac_lutas )) +
geom_text(aes(label = percent(frac_lutas)) ) +
scale_fill_gradient(low = "white", high = "darkgreen", labels = percent ) +
theme_minimal() 
hdi_desc <- read_csv("dados/hdi/desc.csv")
hdi_tudo <- read_csv("dados/hdi/atlas.csv") %>%
mutate(
regiao = case_when(
uf %/% 10 == 1 ~ "N",
uf %/% 10 == 2 ~ "NE",
uf %/% 10 == 3 ~ "SE",
uf %/% 10 == 4 ~ "S",
uf %/% 10 == 5 ~ "CO",
TRUE ~ NA_character_
)
) %>%
select(
espvida,
anos_estudo = e_anosestudo,
gini,
pind,
pren10ricos,
prentrab,
t_banagua,
regiao
)
ggpairs(hdi_tudo) +
theme_minimal()
jogos <- read_csv("dados/football_events/ginf.csv") %>%
mutate(
resultado = case_when(
fthg > ftag ~ "Casa",
fthg < ftag ~ "Fora",
TRUE ~ "Empate"
)
) %>%
count(season, country, resultado)
ggplot(jogos) +
geom_col(aes(y = n, fill = resultado, x = as.factor(season)), position = "fill") +
facet_wrap(~country) +
theme_minimal() +
scale_fill_manual(values = wes_palette("Royal2")) +
scale_y_continuous(labels = percent) +
theme(
axis.text.x = element_text(angle = 90),
legend.position = "top",
text = element_text(family = "Rockwell Condensed")
) +
labs(x = "Temporada", y = "% Vitórias", fill = "Resultado")
jogos <- read_csv("dados/football_events/ginf.csv") %>%
filter(season != 2017) %>%
select(season, country, Fora = ftag, Casa = fthg) %>%
pivot_longer(cols = c(Fora, Casa), names_to = "Time", values_to = "Gols" ) %>%
group_by(season, country, Time) %>%
summarise(Gols = sum(Gols))
ggplot(jogos) +
geom_col(aes(y = Gols, fill = Time, x = as.factor(season)), position = "stack") +
facet_wrap(~country) +
theme_minimal() +
scale_fill_manual(values = wes_palette("Royal2")) +
theme(
axis.text.x = element_text(angle = 90),
legend.position = "top",
text = element_text(family = "Rockwell Condensed")
) +
labs(x = "Temporada", y = "Gols", fill = "Time")
Repare que há muitos países, ou seja, muitas categorias.
Usando a biblioteca forcats e sua função fct_lump conseguimos criar uma categoria de “Outros”
gdps <- wb(indicator = "NY.GDP.MKTP.PP.KD", country = "countries_only")
gdps_tratado <- gdps %>%
mutate(date = as.integer(date) ) %>%
group_by(country) %>%
mutate(ultimo_gdp = last(value, order_by = date)) %>%
ungroup() %>%
mutate(country = fct_lump(country, n = 7, w = ultimo_gdp, other_level = "Outros")) %>%
group_by(country, date) %>%
summarise(value = sum(value),
ultimo_gdp = sum(ultimo_gdp)) %>%
ungroup() %>%
mutate(country = fct_reorder(country, ultimo_gdp ))
ggplot(gdps_tratado ) +
geom_area(aes(x = date, y = value, fill = country), position = "fill") +
theme_minimal() +
scale_fill_brewer(palette = "Accent") +
theme(
axis.text.x = element_text(angle = 90),
legend.position = "top",
text = element_text(family = "CMU Classical Serif")
) +
labs(x = "Ano", y = "PIB PPP", fill = "País")
gdps_tratado <- gdps %>%
mutate(date = as.integer(date) ) %>%
group_by(country) %>%
mutate(ultimo_gdp = last(value, order_by = date)) %>%
ungroup() %>%
mutate(country = fct_lump(country, n = 7, w = ultimo_gdp, other_level = "Outros")) %>%
group_by(country, date) %>%
summarise(value = sum(value),
ultimo_gdp = sum(ultimo_gdp)) %>%
ungroup() %>%
mutate(country = fct_reorder(country, ultimo_gdp ))
ggplot(gdps_tratado ) +
geom_area(aes(x = date, y = value, fill = country) ) +
theme_minimal() +
scale_fill_brewer(palette = "Accent") +
theme(
axis.text.x = element_text(angle = 90),
legend.position = "top",
text = element_text(family = "CMU Serif Extra")
) +
labs(x = "Ano", y = "PIB PPP", fill = "País")
Aqui cabe uma torta se houver poucas categorias
gdps_tratado <- gdps %>%
mutate(date = as.integer(date)) %>%
filter(date == max(date) ) %>%
mutate(country = fct_lump(country, n = 1, w = value, other_level = "Resto")) %>%
group_by(country) %>%
summarise(PIB = sum(value))
ggplot(gdps_tratado, aes(x = "", y = PIB, fill = country)) +
geom_bar( width = 1, stat = "identity", color = "white") +
coord_polar("y", start = 0) +
scale_fill_brewer(palette = "Accent") +
theme_void()
gdps_tratado <- gdps %>%
mutate(date = as.integer(date)) %>%
filter(date == max(date) ) %>%
mutate(country = fct_lump(country, n = 50, w = value, other_level = "Resto")) %>%
group_by(country) %>%
summarise(PIB = sum(value))
ggplot(gdps_tratado, aes( fill = country, area = PIB, label = country)) +
geom_treemap() +
geom_treemap_text(grow = TRUE) +
theme(
legend.position = "none"
)
Pesquisando imagens no Google vemos as percepções a respeito de “Machine Learning” e “Statistical Learning”
(Mello 2019)
Temos acesso a um pedaço dos dados existentes, que usamos para nos ajudar a especificar e treinar o nosso modelo.
Podemos adotar como premissa que existe uma função que traduz como o universo funciona. Esta função determina qual valor \(y\) da variável dependente acontece quando as variáveis independentes assumem valores \(x\), onde \(x\) é um vetor de tamanho \(p\). \(p\) é o número de variáveis dependentes usadas. Esta função não determina perfeitamente \(y\), então temos sempre um \(\epsilon\), um ruído que pode ser devido a variáveis dependentes desconhecidas ou efeitos que não podem ser mensurados.
\[y = f(x) + \epsilon\]
Onde:
\(f(x)\) é uma função desconhecida e \(E(\epsilon) = 0\),
\(\epsilon\) é independente de \(x\) e \(Var(\epsilon) = \sigma^2\) é constante em relação a \(x\).
Nós não temos acesso a esta \(f()\) que representa a VERDADE, mas tentamos estimar uma \(\hat{f}()\)
Nosso modelo de Aprendizado Estatístico é uma fábrica de \(\hat{f}()\) que tenta aproximar a \(f()\) verdadeira.
Existem fábricas para todos os gostos. Desde fábricas que produzem funções simples e inteligíveis até funções do tipo caixa preta.
Esta fábrica de \(\hat{f}()\) recebe como insumo variáveis retiradas do ambiente. Essas variáveis podem ser tratadas de forma a deixar as coisas mais fáceis para o modelo, de modo a deixá-lo na cara do gol como Gérson fazia na copa de 70. Esta etapa de preparar as entradas se chama Feature Engineering. Para fazer isso, usamos tudo que aprendemos sobre manipulação de dados.
A predição é o ato de descobrir \(y\) a partir de \(x\).
A inferência é o entendimento de como valores diferentes de x afetam y
Se rodarmos \(\hat{f}(x_0)\) treinando o modelo com vários conjuntos de treinamento e testasse ele com uma entrada específica qualquer para a qual \(y_0 = f(x_0) + \epsilon\):
\[E[y_o - \hat{f}(x_0)] = Var(\hat{f}(x_0)) + [Bias(\hat{f}(x_0))]^2 + Var(\epsilon)\]
\(Var(\hat{f}(x_0)\) é a variância do modelo: o quanto ele muda quando treinamos ele com outro conjunto de treinamento
\([Bias(\hat{f}(x_0))]^2\) é o erro introduzido por tentarmos aproximar um problema complexo fda vida real com uso de um modelo simplificado.
Os modelos variam muito com relação a “transparência da caixa”
Há uma relação inversamente proporcional (grosso modo) entre:
Flexibilidade do modelo, ou poder de se adaptar a uma relação complexa entre entrada e saída
Interpretabilidade, ou facilidade de o usuário do modelo (piloto) inferir as relações entre as variáveis independentes e as variáveis dependentes
(James et al. 2013)
Além da falta de interpretabilidade, outra questão pode contar contra os modelos mais complexos, com mais poder de aprender nuances sobre os dados: eles podem aprender características específicas dos dados de treinamento.
(Hastie, Tibshirani, and Friedman 2001)
Model Selection: estimação da performance de vários modelos a fim de escolher o melhor
Model Assessment: tendo escolhido o melhor modelo, estimação do erro de generalização em novos dados
Uma forma de usar todos os dados (EXCETO OS DADOS DE TESTE) tanto para treinar quando para validar é revezando que partes dos dados estão sendo usados para treinamento e para validação.
Podemos dividir os dados em \(k\) partes. Cada vez uma parte é escolhida para ser a validação.
(Robot 2019)
Primeiro podemos fazer nossa exploração preliminar.
Como exemplo, podemos executar um pequeno tratamento das features usando a função pivot_longer() de um jeito diferente, só possível na últiam versão
ufc_data <- read_csv("dados/ufc/data.csv")
ufc_raw_data <- read_csv2("dados/ufc/raw_total_fight_data.csv")
ufc_raw_fighter <- read_csv2("dados/ufc/raw_fighter_details.csv")
ufc_raw_pre <- read_csv2("dados/ufc/preprocessed_data.csv")
ufc_cada_lutador <- bind_cols(ufc_data, ufc_raw_data) %>%
pivot_longer(
cols = matches("[BR]_.*"),
names_to = c("lutador",".value") ,
names_sep = "_"
) %>%
mutate(
ganhou = str_sub(Winner,1,1) == lutador
) %>%
mutate(aprov = wins/(wins+losses) )
ggplot(ufc_cada_lutador) +
geom_boxplot(aes(x = ganhou, y = aprov )) +
facet_wrap(~weight_class) +
theme_minimal()
ufc_aprend <- bind_cols(ufc_data, ufc_raw_data) %>%
filter(Winner != "Draw") %>%
mutate(
R_aprov = R_wins / (R_wins + R_losses),
B_aprov = B_wins / (B_wins + B_losses),
) %>%
select(
Winner,
weight_class,
B_Weight_lbs,
B_age,
R_Weight_lbs,
R_age,
B_age,
R_aprov,
B_aprov
) %>%
mutate(
winner_color = Winner,
zebra_blue = if_else(Winner == "Blue",1,0),
red_mais_velho = R_age - B_age,
idade_red = R_age
) %>%
filter(
!is.na(idade_red) & !is.na(red_mais_velho) & !is.na(zebra_blue)
)Como exemplo, podemos observar a relação entre a diferença de idade dos lutadores e o vencedor da luta.
Normalmente o lutador vermelho é o favorito. mas podemos ver que quando o lutador vermelho é velho e mais velho que o azul, ele costuma perder.
ggplot(ufc_aprend ) +
geom_jitter(aes( y = red_mais_velho, x = R_age, color = winner_color ), alpha = 0.15) +
scale_color_manual( values =c("Blue" = "blue", "Red" = "red")) +
facet_wrap( ~weight_class ) +
theme_minimal()
purrPreparando um grid das variáveis independentes
red_mais_velho <- seq(
from = min(ufc_aprend$red_mais_velho, na.rm = TRUE),
to = max(ufc_aprend$red_mais_velho, na.rm = TRUE),
length.out = 50
) %>%
enframe(name = "1", value = "red_mais_velho")
idade_red <- seq(
from = min(ufc_aprend$idade_red, na.rm = TRUE),
to = max(ufc_aprend$idade_red, na.rm = TRUE),
length.out = 50
) %>%
enframe(name = "2", value = "idade_red")
weight_classes <- ufc_aprend %>%
select(weight_class) %>%
distinct()
grid_previsao <- red_mais_velho %>%
crossing(idade_red) %>%
crossing(weight_classes)
ggplot(grid_previsao %>% filter(weight_class == "Heavyweight")) +
geom_point(aes(x = idade_red, y = red_mais_velho ), size = 0.01) +
theme_minimal()
purrRodando o modelo com o grid
previsoes <- augment(modelo_lm, newdata = grid_previsao )
#idade_red + red_mais_velho
ggplot(ufc_aprend ) +
geom_jitter(aes( y = red_mais_velho, x = idade_red, color = winner_color ), alpha = 0.3) +
scale_color_manual( values =c("Blue" = "blue", "Red" = "red")) +
facet_wrap( ~weight_class ) +
geom_point(
data = previsoes %>% filter(.fitted > 0.5),
aes(y = red_mais_velho, x = idade_red ),
color = "lightblue",
alpha = 0.1
) +
geom_point(
data = previsoes %>% filter(.fitted < 0.5),
aes(y = red_mais_velho, x = idade_red ),
color = "lightcoral",
alpha = 0.1
) +
theme_minimal()
grid_previsao_com_modelo <- grid_previsao %>%
group_by(weight_class) %>%
nest() %>%
inner_join(modelo_ufc_por_classe, by = c("weight_class")) %>%
mutate(previsoes = map2(.x = modelo, .y = data, .f = ~augment(x = .x, newdata = .y) )) %>%
unnest_legacy(previsoes)
ggplot(ufc_aprend ) +
geom_jitter(aes( y = red_mais_velho, x = idade_red, color = winner_color ), alpha = 0.3) +
scale_color_manual( values =c("Blue" = "blue", "Red" = "red")) +
facet_wrap( ~weight_class ) +
geom_point(
data = grid_previsao_com_modelo %>% filter(.fitted > 0.5),
aes(y = red_mais_velho, x = idade_red ),
color = "lightblue",
alpha = 0.1
) +
geom_point(
data = grid_previsao_com_modelo %>% filter(.fitted < 0.5),
aes(y = red_mais_velho, x = idade_red ),
color = "lightcoral",
alpha = 0.1
) +
theme_minimal()
#ALGO ERRADO AO RODAR O MARKDOWN. NÃO SEI
# modelo_ufc_por_classe_treino_gam <- ufc_aprend_treino %>%
# # group_by(weight_class) %>%
# nest(data =c(zebra_blue, idade_red, red_mais_velho ) ) %>%
# mutate(
# modelo = map(data, ~gam(formula = zebra_blue ~ s(idade_red) + s(red_mais_velho), data = .x))
# ) %>%
# select(-data)
#
#
# resultados_teste_gam <- ufc_aprend_teste %>%
# group_by(weight_class) %>%
# nest_legacy() %>%
# inner_join(modelo_ufc_por_classe_treino_gam, by = c("weight_class") ) %>%
# mutate(
# previsao = map2(.x = modelo, .y = data, .f = ~predict(.x, .y) )
# ) %>%
# unnest(cols = c("data","previsao")) %>%
# ungroup() %>%
# mutate(
# actual = if_else(zebra_blue == 1, "B", "R"),
# predicted = if_else(previsao > 0.5, "B", "R"),
# ) %>%
# mutate(
# actual = fct_relevel(actual, "R", "B"),
# predicted = fct_relevel(predicted, "R", "B")
# ) %>%
# select(
# actual,
# predicted
# )
#
# write_rds(resultados_teste_gam, "dados/ufc/resultados_teste_gam.rds")
resultados_teste_gam <- read_rds("dados/ufc/resultados_teste_gam.rds")
confusionMatrix(resultados_teste_gam$predicted, resultados_teste_gam$actual )## Confusion Matrix and Statistics
##
## Reference
## Prediction R B
## R 716 310
## B 48 53
##
## Accuracy : 0.6823
## 95% CI : (0.6543, 0.7095)
## No Information Rate : 0.6779
## P-Value [Acc > NIR] : 0.3884
##
## Kappa : 0.1026
##
## Mcnemar's Test P-Value : <2e-16
##
## Sensitivity : 0.9372
## Specificity : 0.1460
## Pos Pred Value : 0.6979
## Neg Pred Value : 0.5248
## Prevalence : 0.6779
## Detection Rate : 0.6353
## Detection Prevalence : 0.9104
## Balanced Accuracy : 0.5416
##
## 'Positive' Class : R
## Plotando as regiões deste modelo mais complexo
#
# grid_previsao_com_modelo_gam <- grid_previsao %>%
# group_by(weight_class) %>%
# nest() %>%
# inner_join(modelo_ufc_por_classe_treino_gam, by = c("weight_class")) %>%
# mutate(previsoes = map2(.x = modelo, .y = data, .f = ~predict( .x, .y) )) %>%
# select(-modelo) %>%
# unnest(cols = c("previsoes","data")) %>%
# ungroup()
#
#
# write_rds(grid_previsao_com_modelo_gam, "dados/ufc/grid_previsao_com_modelo_gam.RDS")
grid_previsao_com_modelo_gam <- read_rds("dados/ufc/grid_previsao_com_modelo_gam.RDS")
ggplot( bind_rows(ufc_aprend_classes_selec) ) +
geom_jitter(aes( y = red_mais_velho, x = idade_red, color = winner_color ), alpha = 0.3) +
scale_color_manual( values =c("Blue" = "blue", "Red" = "red")) +
facet_wrap( ~weight_class ) +
geom_point(
data = grid_previsao_com_modelo_gam %>% filter(previsoes > 0.5),
aes(y = red_mais_velho, x = idade_red ),
color = "lightblue",
alpha = 0.1
) +
geom_point(
data = grid_previsao_com_modelo_gam %>% filter(previsoes < 0.5),
aes(y = red_mais_velho, x = idade_red ),
color = "lightcoral",
alpha = 0.1
) +
theme_minimal()
A curva ROC (Receiver Operating Characteristics) foi inventada na época da Segunda Guerra Mundial para avaliar se os operadores de radar americanos estavam detectando confiavelmente aeronaves japonesas a partir de sinais de radar.
A curva mostra, para vários thresholds, qual a fração de verdadeiros positivos (ou sensibilidade) e a fração de falsos positivos (fall-out, ou \(1 - especificidade\) ).
Uma métrica numérica que traduz o quão bom um modelo de classificação é consiste na área embaixo desta curva (AUC, Area Under the Curve). Note que quanto mais perto de um essa área, menor a taxa de falsos positivos e maior a sensibilidade
ggplot(model_logistic$pred, aes(m = M, d = obs, color = Resample )) +
geom_roc( labels = FALSE ) +
coord_equal() + style_roc() + ggtitle("ROC", subtitle = "Métricas para diversos thresholds" )+ style_roc()
Cada conjunto de hiperparâmetros do modelo rodado nesta possui dois valores. Um dos valores, \(\alpha\) define que peso deve ser dado a cada tipo de penalização Ridge ou Lasso, onde \(\alpha = 0\) é Ridge puro e \(\alpha = 1\) é Lasso puro. O outro parâmetro, \(\lambda\), regula a intensidade da penalização.
Os resultados para várias parametrizações é mostrado abaixo.
Veja como podemos passar um grid de hiperparâmetros e avaliar esses resultados.
set.seed(13)
model_net <- train(
Diagnosis ~ .,
treino,
metric = "ROC",
method = "glmnet",
trControl = controle_cv,
preProcess = c( "center", "scale"),
tuneGrid = expand.grid(
alpha = c(0,0.5,0.75,0.1,0.125,0.15,0.25,0.5, 1),
lambda = 0:10/500
)
)
plot(model_net)
Os resultados de cada parametrização são mostrados em ordem decrescente de AUC.
resultados_net <- model_net$resample %>%
group_by(alpha, lambda) %>%
summarise(media = mean(ROC), "Desvio-padrão AUC" = sd(ROC)) %>%
arrange(desc(media)) %>%
rename("Média AUC" = media)
head(resultados_net, 20)## # A tibble: 20 x 4
## # Groups: alpha [6]
## alpha lambda `Média AUC` `Desvio-padrão AUC`
## <dbl> <dbl> <dbl> <dbl>
## 1 1 0.02 0.994 0.00746
## 2 1 0.014 0.994 0.00762
## 3 1 0.018 0.994 0.00756
## 4 1 0.016 0.994 0.00761
## 5 1 0.012 0.994 0.00776
## 6 1 0.01 0.994 0.00782
## 7 1 0.008 0.994 0.00777
## 8 0.75 0.006 0.994 0.00790
## 9 0.1 0.01 0.994 0.00715
## 10 0.25 0.01 0.994 0.00732
## 11 0.125 0.018 0.994 0.00732
## 12 0.15 0.008 0.994 0.00720
## 13 0.1 0.02 0.994 0.00730
## 14 1 0.006 0.994 0.00784
## 15 0.1 0.016 0.994 0.00724
## 16 0.1 0.018 0.994 0.00726
## 17 0.15 0.014 0.994 0.00728
## 18 0.1 0.012 0.994 0.00716
## 19 0.75 0.004 0.994 0.00766
## 20 0.125 0.008 0.994 0.00721Uma outra forma de reduzir a dimensionalidade do probelam é usar o método PCA para transformar as variáveis em componentes principais em menor número mas com boa parte da informação original.
veja como pedimos esse pré-processamento.
set.seed(13)
model_net_pca <- train(
Diagnosis ~ .,
treino,
metric = "ROC",
method = "glmnet",
trControl = controle_cv,
preProcess = c( "center", "scale", "pca"),
tuneGrid = expand.grid(
alpha = c(0,0.5,0.75,0.1,0.125,0.15,0.25,0.5, 1),
lambda = 0:10/500
)
)
plot(model_net_pca)
resultados_net_pca <- model_net_pca$resample %>%
group_by(alpha, lambda) %>%
summarise(media = mean(ROC), "Desvio-padrão AUC" = sd(ROC)) %>%
arrange(desc(media)) %>%
rename("Média AUC" = media)
resultados_net_pca## # A tibble: 88 x 4
## # Groups: alpha [8]
## alpha lambda `Média AUC` `Desvio-padrão AUC`
## <dbl> <dbl> <dbl> <dbl>
## 1 1 0.01 0.993 0.00664
## 2 1 0.012 0.993 0.00658
## 3 1 0.014 0.993 0.00648
## 4 1 0.008 0.993 0.00692
## 5 0.75 0.016 0.993 0.00709
## 6 0.75 0.012 0.993 0.00734
## 7 0.75 0.014 0.993 0.00719
## 8 0.75 0.006 0.993 0.00754
## 9 0.75 0.008 0.993 0.00758
## 10 1 0.006 0.993 0.00703
## # ... with 78 more rowsA árvore de decisão particiona o espaço formado pelas variáveis explicativas em subespaços baseando-se na “pureza” desses subespaços com relação à variável dependente.
model_tree <- train(
Diagnosis ~ .,
treino,
metric = "ROC",
method = "rpart",
trControl = controle_cv,
preProcess = c( "center", "scale")
)
resultados_tree <- model_tree$resample %>%
group_by(cp) %>%
summarise(media = mean(ROC), "Desvio-padrão AUC" = sd(ROC)) %>%
arrange(desc(media)) %>%
rename("Média AUC" = media)
resultados_tree## # A tibble: 3 x 3
## cp `Média AUC` `Desvio-padrão AUC`
## <dbl> <dbl> <dbl>
## 1 0.0214 0.942 0.0273
## 2 0.0500 0.932 0.0293
## 3 0.829 0.773 0.193A forma como a árvore de decisão é criada faz com que ela tenha muita variância.
Cada decisão de particionamento é tomada a partir de características que podem ser muito específicas aos dados de treinamento.
Uma ideia usada nas Random Forests é criar conjuntos de treinamento criados a partir do conjunto original, mas retirando amostras de mesmo tamanho COM REPOSIÇÃO. Além disso, cada vez que a árvore é particionada, a partição só pdoe acontecer em \(m\) das variáveis explicativas.
O resultado final é uma média do resultado dessas várias árvores
Estas duas mudanças fazem com que o modelo tenha uma variância muito menor do que as árvore de decisão simples
model_ranger <- train(
Diagnosis ~ .,
treino,
metric = "ROC",
method = "ranger",
trControl = controle_cv,
preProcess = c( "center", "scale"),
verbose = TRUE,
tuneGrid = expand.grid(
mtry = seq(2, by = 2, to = 20),
splitrule = c("gini", "extratrees"),
min.node.size = 1
)
)
plot(model_ranger)
resultados_ranger <- model_ranger$resample %>%
group_by(mtry, splitrule) %>%
summarise(media = mean(ROC), "Desvio-padrão AUC" = sd(ROC)) %>%
arrange(desc(media)) %>%
rename("Média AUC" = media)
resultados_ranger## # A tibble: 20 x 4
## # Groups: mtry [10]
## mtry splitrule `Média AUC` `Desvio-padrão AUC`
## <dbl> <fct> <dbl> <dbl>
## 1 2 extratrees 0.992 0.00606
## 2 4 extratrees 0.992 0.00647
## 3 8 extratrees 0.992 0.00682
## 4 6 extratrees 0.991 0.00686
## 5 16 extratrees 0.991 0.00778
## 6 10 extratrees 0.991 0.00759
## 7 12 extratrees 0.990 0.00808
## 8 14 extratrees 0.990 0.00847
## 9 20 extratrees 0.990 0.00848
## 10 18 extratrees 0.990 0.00913
## 11 2 gini 0.989 0.00933
## 12 6 gini 0.988 0.0104
## 13 18 gini 0.988 0.0101
## 14 16 gini 0.988 0.0101
## 15 10 gini 0.988 0.0104
## 16 4 gini 0.988 0.0105
## 17 20 gini 0.988 0.0101
## 18 8 gini 0.987 0.0106
## 19 14 gini 0.987 0.0104
## 20 12 gini 0.987 0.0108É possível propor uma forma de avaliação do impacto das variáveis indireta e agnóstica ao modelo, ou seja, que pode ser usada da mesma forma para qualquer modelo e .
Para avaliar o impacto de cada variável nos resultados de cada um dos 5 modelos, podemos gerar casos sintéticos e avaliar a probabilidade de esses casos serem malignos segundo cada modelo.
Dois grupos de casos sintéticos foram gerados.
No primeiro grupo, para cada variável explicativa, foram criados casos onde o valor das outras variáveis.
No segundo grupo, para avaliar uma região onde os casos têm mais chance de ser malignos, as variáveis que não estão sendo avaliadas a cada gráfico são mantidas no quantil 0.65.
medias_treino <- treino %>%
select(-Diagnosis) %>%
summarise_all(mean) %>%
gather(variavel, media)
dps_treino <- treino %>%
select(-Diagnosis) %>%
summarise_all(sd) %>%
gather(variavel, dp)
mediana_treino <- treino %>%
select(-Diagnosis) %>%
summarise_all(median) %>%
gather(variavel, mediana)
variaveis <- medias_treino %>%
select(variavel)
caso_media <- medias_treino %>%
spread(variavel, media)
quantis <- treino %>%
select(-Diagnosis) %>%
gather(variavel, valor) %>%
crossing( tibble( q = seq(0,1, 0.05)) ) %>%
group_by(variavel,q ) %>%
summarise( valor = quantile(valor, mean(q)))
casos <- variaveis %>%
crossing(variaveis %>% rename(fixa = variavel)) %>%
crossing(tibble( q = seq(0,1, 0.05)) ) %>%
arrange(q, variavel, fixa) %>%
mutate(caso = (row_number()-1) %/% nrow(variaveis) + 1) %>%
mutate(q_cada_variavel = if_else(variavel == fixa, q, 0.5 )) %>%
inner_join(quantis, by=c("q_cada_variavel" = "q", "fixa" = "variavel") ) %>%
select(-q_cada_variavel) %>%
spread(fixa, valor)
prev_log <- predict(model_logistic, casos, type = "prob") %>%
mutate(modelo = "Logistica") %>%
bind_cols(casos)
prev_ranger <- predict(model_ranger, casos, type = "prob") %>%
mutate(modelo = "Floresta") %>%
bind_cols(casos)
prev_pca <- predict(model_net_pca, casos, type = "prob") %>%
mutate(modelo = "Lasso-Ridge PCA") %>%
bind_cols(casos)
prev <-
bind_rows(prev_log, prev_ranger, prev_pca)
prev %>%
ggplot() +
geom_line( aes(x = q, y = M, color = modelo) ) +
facet_wrap( ~variavel) +
theme_minimal() +
theme(
aspect.ratio = .4,
strip.text = element_text(size = 5),
strip.background.y = element_rect(size = c(0.1,0.1)),
strip.switch.pad.grid = unit(.01,"cm"),
strip.switch.pad.wrap = unit(.01,"cm"),
axis.text = element_text(size = 4)
) +
theme(legend.position = "top" )
casos <- variaveis %>%
crossing(variaveis %>% rename(fixa = variavel)) %>%
crossing(tibble( q = seq(0,1, 0.05)) ) %>%
arrange(q, variavel, fixa) %>%
mutate(caso = (row_number()-1) %/% nrow(variaveis) + 1) %>%
mutate(q_cada_variavel = if_else(variavel == fixa, q, 0.65 )) %>%
inner_join(quantis, by=c("q_cada_variavel" = "q", "fixa" = "variavel") ) %>%
select(-q_cada_variavel) %>%
spread(fixa, valor)
prev_log <- predict(model_logistic, casos, type = "prob") %>%
mutate(modelo = "Logistica") %>%
bind_cols(casos)
prev_ranger <- predict(model_ranger, casos, type = "prob") %>%
mutate(modelo = "Floresta") %>%
bind_cols(casos)
prev_pca <- predict(model_net_pca, casos, type = "prob") %>%
mutate(modelo = "Lasso-Ridge PCA") %>%
bind_cols(casos)
prev <-
bind_rows(prev_log, prev_ranger, prev_pca)
prev %>%
ggplot() +
geom_line( aes(x = q, y = M, color = modelo) ) +
facet_wrap( ~variavel) +
theme_minimal() +
theme(
aspect.ratio = .4,
strip.text = element_text(size = 5),
strip.background.y = element_rect(size = c(0.1,0.1)),
strip.switch.pad.grid = unit(.01,"cm"),
strip.switch.pad.wrap = unit(.01,"cm"),
axis.text = element_text(size = 4)
) +
theme(legend.position = "top" )
roc_best_modelo <- function(modelo)
{
dados_roc <- modelo$pred %>%
semi_join(modelo$bestTune) %>%
mutate(obs = if_else(obs == "B", 0, 1))
dados_curva <- roc(dados_roc$obs, dados_roc$M)
tibble( sensitivities = dados_curva$sensitivities,
specificities= dados_curva$specificities,
thresholds = dados_curva$thresholds,
auc = dados_curva$auc
)
}
pred_best_modelo <- function(modelo)
{
modelo$pred %>%
semi_join(modelo$bestTune)
}
plota_roc_modelo <- function(modelo, prob, nome)
{
roc_modelo <- roc_best_modelo(modelo)
pred_modelo <- pred_best_modelo(modelo)
threshold_escolhido <- roc_modelo %>%
filter(thresholds >= prob) %>%
top_n(n = 1, wt = desc(thresholds ))
ggplot(pred_modelo ) +
geom_roc( labels = FALSE, labelround = 2, labelsize = 3, n.cuts = 0, aes(m = M, d = obs ) ) +
geom_point(data = threshold_escolhido, aes(y = sensitivities, x = 1-specificities), size = 3) +
geom_text(data = threshold_escolhido, aes(y = sensitivities, x = 1-specificities, label = paste0(percent(sensitivities), " / ", percent(1-specificities)) ), nudge_x = .11, nudge_y = -.05, size = 3 ) +
coord_equal() + style_roc() + ggtitle(paste0("ROC - ", nome), subtitle = paste0("Corte na probabilidade ", prob )) + style_roc()
}
plota_roc_ensemble <- function(modelo, prob, nome)
{
prob <- 0.3
nome <- "Ensemble"
dados_curva <- roc(pred_ensemble_mediana$obs, pred_ensemble_mediana$M)
roc_modelo <- tibble( sensitivities = dados_curva$sensitivities,
specificities= dados_curva$specificities,
thresholds = dados_curva$thresholds,
auc = dados_curva$auc
)
pred_modelo <- pred_ensemble_mediana
threshold_escolhido <- roc_modelo %>%
filter(thresholds >= prob) %>%
top_n(n = 1, wt = desc(thresholds ))
ggplot(pred_modelo ) +
geom_roc( labels = FALSE, labelround = 2, labelsize = 3, n.cuts = 0, aes(m = M, d = obs ) ) +
geom_point(data = threshold_escolhido, aes(y = sensitivities, x = 1-specificities), size = 3) +
geom_text(data = threshold_escolhido, aes(y = sensitivities, x = 1-specificities, label = paste0(percent(sensitivities), " / ", percent(1-specificities)) ), nudge_x = .11, nudge_y = -.05, size = 3 ) +
coord_equal() + style_roc() + ggtitle(paste0("ROC - ", nome), subtitle = paste0("Corte na probabilidade ", prob )) + style_roc()
}
plota_roc_modelo(model_logistic, 1e-8, "Logistica")
plota_roc_ensemble <- function(modelo, prob, nome)
{
prob <- 0.17
nome <- "Ensemble"
dados_curva <- roc(pred_ensemble_mediana$obs, pred_ensemble_mediana$M)
roc_modelo <- tibble( sensitivities = dados_curva$sensitivities,
specificities= dados_curva$specificities,
thresholds = dados_curva$thresholds,
auc = dados_curva$auc
)
pred_modelo <- pred_ensemble_mediana
threshold_escolhido <- roc_modelo %>%
filter(thresholds >= prob) %>%
top_n(n = 1, wt = desc(thresholds ))
ggplot(pred_modelo ) +
geom_roc( labels = FALSE, labelround = 2, labelsize = 3, n.cuts = 0, aes(m = M, d = obs ) ) +
geom_point(data = threshold_escolhido, aes(y = sensitivities, x = 1-specificities), size = 3) +
geom_text(data = threshold_escolhido, aes(y = sensitivities, x = 1-specificities, label = paste0(percent(sensitivities), " / ", percent(1-specificities)) ), nudge_x = .11, nudge_y = -.05, size = 3 ) +
coord_equal() + style_roc() + ggtitle(paste0("ROC - ", nome), subtitle = paste0("Corte na probabilidade ", prob )) + style_roc()
}
plota_roc_ensemble()
dados_confusao_teste_todos_tidy <- dados_confusao_teste_todos %>%
gather(atributo, valor, - pred,-Diagnosis, - modelo, - linha)
dados_confusao_teste_todos_certos <-
dados_confusao_teste_todos_tidy %>%
filter(pred == Diagnosis)
dados_confusao_teste_todos_falso_negativo <-
dados_confusao_teste_todos_tidy %>%
filter(pred == "M" & Diagnosis == "B")
dados_confusao_teste_todos_falso_positivo <-
dados_confusao_teste_todos_tidy %>%
filter(pred == "B" & Diagnosis == "M")
ggplot() +
geom_jitter(data = dados_confusao_teste_todos_certos %>% filter(Diagnosis == "M"), aes(x = valor, y = 0), alpha = 0.05, size = 0.5, color = "red") +
geom_jitter(data = dados_confusao_teste_todos_certos %>% filter(Diagnosis == "B") , aes(x = valor, y = 0), alpha = 0.05, size = 0.5, color = "blue") +
geom_jitter(data = dados_confusao_teste_todos_falso_negativo, aes(x = valor, y = 0), alpha = 0.5, color = "blue", size = 0.5) +
geom_jitter(data = dados_confusao_teste_todos_falso_positivo , aes(x = valor, y = 0), alpha = 0.5, color = "red", size = 0.5) +
facet_wrap(~atributo, scales = "free") +
theme_minimal() +
theme(
aspect.ratio = .4,
strip.text = element_text(size = 5),
strip.background.y = element_rect(size = c(0.1,0.1)),
strip.switch.pad.grid = unit(.01,"cm"),
strip.switch.pad.wrap = unit(.01,"cm"),
axis.text = element_text(size = 4) ,
axis.text.y = element_blank()
)